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waves.  This  report  refines  and  extends  a  numerical  method  based  on  Fourier  series  solutions  of  the 
governing  equations.  The  resultant  solitary  wave  properties  agree  with  those  found  by  independent 
methods  to  several  significant  figures,  except  for  waves  having  amplitudes  within  a  few  percent  of 
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that  of  the  highest  wave.  There  appears  to  be  little  degradation  in  the  method’s  accuracy  for 
the  highest  and  very  high  waves,  when  compared  to  waves  of  intermediate  height.  This  report 
discusses  possible  reasons  for  the  failure  of  other  numerical  work  to  achieve  high  accuracy. 
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HIGH  SOLITARY  WAVES  IN  WATER: 
A  REFINED  NUMERICAL  METHOD 


1.  INTRODUCTION 

This  report  describes  a  numerical  method  for  analyzing  the  properties  of  solitary  waves  in  water 
and  assesses  the  method's  accuracy.  A  second  report  (J.  M.  Witting,  "High  Solitary  Waves  in  Water: 
Results  of  Calculations,"  NRL  Report  8505)  presents  detailed  and  comprehensive  results. 

Figure  1  presents  the  motivation  for  the  research  and  hints  at  some  of  the  results.  The  figure 
plots  the  square  of  the  speed  of  a  solitary  wave  in  water  against  its  amplitude.  The  plot  is  nondimen- 
sional;  speeds  are  referenced  to  the  speed  of  long  linear  waves  in  water;  amplitudes  are  referenced  to 
the  water  depth  well  away  from  the  solitary  wave.  To  make  the  figure  legible,  only  high  solitary  waves 
are  shown.  The  complete  curve  would  extend  downward  to  a  speed  of  unity  and  an  amplitude  of  zero. 
It  terminates  with  the  highest  wave  in  water.  The  highest  wave  has  a  "kineticity”  twice  the  amplitude.* 
Figure  1  shows  the  limiting  relationship  as  a  barrier. 


azti 


Fig.  I  —  Amplitude-speed  relationship  Displayed  are  sample  data  from  several  theories  that  claim 
high  accuracy.  The  numerical  results  of  previous  investigators  generally  fall  below  those  obtained 
by  other  means,  except  for  the  highest  wave.  Our  calculations  of  the  relationship  data  agree  with 
those  of  Longuet-Higgins  and  Fenton  (1974)  and  with  Byatt-Smith  and  Longuet-HIggins  (1976), 
except  for  very  high  waves. 


Manuscript  submitted  on  June  24,  1981 

'This  fact  has  been  known  since  Stokes  (see  Lamb.  1932  Art  250),  if  not  before,  it  is  the  result  of  the  highest  wave's  possessing 
a  sharp  crest,  i.e.,  a  discontinuous  slope  The  only  way  that  this  can  occur  is  if  the  fluid  at  the  crest  moves  at  the  same  speed  as 
the  wave  itself  The  Bernoulli  Law  then  demands  the  limiting-height-speed  relationship. 
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Methods  used  to  deduce  the  properties  of  solitary  waves  fall  into  three  classes:  (a)  expansions 
about  a  small  parameter,  e  g.  the  amplitude,  (b)  numerical  methods,  and  (c)  others.  No  matter  what 
the  method,  reliance  on  numerical  computation  is  necessary  to  approach  the  higher  waves.  Figure  I 
presents  the  results  of  our  calculation,  and  displays  sample  results  from  each  of  the  three  methods 
whose  authors  claim  4  to  5  figure  accuracy.  Prior  to  this  work,  mutual  agreement  of  the  stated  accura¬ 
cies  could  be  found  only  between  the  results  of  Longuet-Higgins  and  Fenton  (1974)  and  those  of 
Byalt-Smith  and  Longuet-Higgins  (1976).  Except  for  the  highest  wave,  the  numerical  methods  pro¬ 
duced  somewhat  smaller  speeds  for  a  given  amplitude  than  did  the  others.  All  previous  works  claiming 
high  accuracy  give  a  limiting  wave  height  within  0.0015  of  0.8270. 

Figure  1  shows  that  for  all  but  the  very  highest  waves,  perfect  (to  the  resolution  of  the  figure) 
agreement  is  found  between  the  numerical  method  presented  here  (the  line),  and  the  results  of 
Longuet-Higgins  and  Fenton  (1974)  and  of  Byatt-Smith  and  Longuet-Higgins  (1976).  The  dilemma 
remains,  however,  that  the  method  described  here  also  yields  a  highest  wave  higher  than  previously 
reported.  We  feel  that  a  detailed  description  of  our  numerical  method  is  called  for  to  discuss  why  our 
method  gives  more  accurate  results  than  other  numerical  methods  (Miles,  1980  reviews  the  status  of 
solitary  wave  research  and  points  out  the  unsatisfactory  status  of  knowledge  of  high  solitary  waves). 

This  report  is  organized  as  follows:  Section  2  develops  the  basic  theory  for  describing  a  solitary 
wave  and  defines  certain  conformal  transformations  that  are  useful  for  numerical  calculations.  Section 
3  describes  the  representation  of  the  solution  and  identifies  features  of  the  solitary  wave  which  are  com¬ 
puted.  Section  4  describes  the  numerical  procedures  used  in  obtaining  solutions.  Section  5  compares 
the  results  to  those  of  earlier  investigators.  Section  6  briefly  summarizes  the  results  and  offers  conclu¬ 
sions. 

2.  BASIC  THEORY  AND  USEFUL  TRANSFORMATIONS 

A.  The  Kind  of  Solitary  Wave  Considered 

Figure  2a  shows  an  idealization  of  solitary  waves  in  water.  Apart  from  dissipative  effects,  the 
wave  propagates  at  a  constant  speed  and  without  change  of  form  over  an  impermeable  and  horizontal 
bottom.  The  conventional  mathematical  model  describing  the  solitary  wave  takes  the  fluid  to  be 
incompressible  and  inviscid,  with  all  fluid  motions  irrotational.  In  addition,  the  conventional  model 
neglects  surface  stresses  imposed  by  the  relatively  undense  air,  and  imposes  a  constant  pressure  boun¬ 
dary  condition  at  the  air-water  interface.  Because  only  pressure  gradients  appear  in  the  equations  of 
motion,  the  upper  surface  of  the  water  is  considered  stress-free.  The  model  contemplates  only  plane 
waves,  so  that  motion  is  confined  to  the  x-y  plane  (see  Fig.  2a).  By  assuming  at  the  start  that  the  wave 
moves  at  constant  speed  and  shape  through  still  water,  one  may  make  a  Galilean  transformation  to 
wave  coordinates  where  the  wave  is  stationary.  Viewed  from  wave  coordinates  the  situation  is  then  one 
of  a  steady,  incompressible,  irrotational  flow  in  two  dimensions. 

Two  equations  completely  govern  the  local  motion  in  the  interior  of  the  fluid.  These  are  that 
div  u  -  0  (incompressibility)  and  (curlu),  —  0  (irrotationality).  It  follows  that  a  velocity  potential  <t> 
and  a  stream  function  are  tied  to  Cauchy-Riemann  conditions'. 

u  -  b<t>/bx  -  9t It/dy, 

v  -  b<f>/by  —  -di Ir/bx. 

As  is  customary,  define  a  pair  of  complex  variables  z  =  x  +  iy  and  tv  —  <b  +  ii/i.  Then  the  local  motion 
in  the  interior  of  the  fluid  is  satisfied  by  any  analytic  function  w  —  w(z),  or,  alternately  z  =  z(»).  Fig¬ 
ure  2b  is  the  map  to  tvspace. 

It  remains  to  pose  boundary  conditions.  These  are  the  following:  the  flow  at  the  bottom  of  the 
channel  is  horizontal;  flows  away  from  the  wave  are  horizontal  and  shear-free,  and  the  same  on  each 
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side  of  the  wave;  the  pressure  along  the  free  surface  vanishes.  Identifying  this  upper  boundary  in 
advance  is  possible  in  w-space  (it  corresponds  to  1 b  =  constant).  It  is  not  possible  in  z-space.  There¬ 
fore,  most  researchers  since  Stokes  (1880)  have  selected  w  to  be  the  independent  variable,  then 
z  =  z(m  )  is  the  desired  solution.  With  this  choice  of  independent  variable  the  boundary  conditions  are. 


0 

II 

:x 

along  tli  =  0 

(2) 

dz  _1_ 

dw  U 

as  <£  —  ±  00 

(3) 

gy  +  ^n-^-gh  +  ~u7 

2  dw  2 

along  <1/  =  ib0 

(4) 

where  g  is  the  acceleration  of  gravity,  h  is  the  water  depth  away  from  the  wave,  and  U  is  the  speed  of 
the  wave  relative  to  still  water.  Like  Yamada  (1957,  1958),  we  nondimensionalize  in  the  way  that 
minimizes  external  parameters  appearing  in  z  and  w.  This  is  the  system  h  —  U  -  1.  Equations  (2), 
(3),  and  (4)  then  become: 


y  -  0  along  <b  =  0  (5) 

— — •  I  as  — *  ±  °°  (6) 

dw 

v-  +  \rF2l  l-^-l3  =  1  +  -i-T2  along  \b  «  I  (7) 

2  dw  2 

where  F  S  U/Jgh  is  the  Froude  Number. 

Henceforth  in  this  report,  solitary  waves  are  defined  to  be  z  -  z(w)  subject  to  Eqs.  (5),  (6),  and 
(7).  As  Fig  I  shows,  there  may  be  more  than  one  solitary  wave  that  corresponds  to  a  given  F.  This 
remarkable  discovery  was  made  by  Longuet-Higgins  and  Fenton  (1974),  and  has  been  confirmed  by 
other  investigators. 

B.  A  Transformation  of  the  Free  Surface  to  a  Unit  Circle 

For  the  purposes  of  performing  a  trustworthy  numerical  analysis,  we  found  it  desirable  to  move 
the  domain  of  the  nasty  boundary  condition,  Eq.  (7),  to  the  unit  circle  by  a  conformal  transformation. 
Levi-Civita  (1925)  and  Yamada  (1957,  1958)  show  how  to  do  this.  The  transformation 

C,  =  -tanh2  w  (8) 

maps  the  region  occupied  by  the  fluid  into  the  interior  of  the  unit  circle  (see  Figs,  lb,  lc). 


The  flow  can  be  described  with  a  new  dependent  variable: 


J2  =  /  In (dtv/dz) 

=  9  +  ir. 

(9) 

In  Eq.  (9)  H  is  the  angle  of  surface  inclination 

and  r  =  In  </,  where  q  is  the  fluid  speed 

dw 

dz 

Referring  to  Fig.  2c,  the  boundary  conditions  are. 

0  -  0 

along  0,  ±  / 

(10) 

a  -  0 

at  ±  / 

(11) 

2  da  1  1 

ql  — cos—  <t  —  — -r  sinS 

d<r  2  ttF‘ 

on  the  unit  circle. 

(12) 
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The  last  condition  is  not  self-evident  (Yamada  1958  gives  a  derivation).  Its  integral  forms  an 
alternative  to  Eq.  (12): 

J  ptr  | 

<j3(<t)  =  o3(0)  H - r  I  sin  0  sec  —<t d<r  on  the  unit  circle.  (13) 

itF2  Jo  2  ; 

! 

C.  Another  Useful  Transformation 

The  transformation  described  in  Section  2B  of  this  report  has  the  unfortunate  side-effect  in  which 
the  point  marked  ±  /  in  Fig.  2c  is  singular.  Although  we  shall  use  numerical  methods  that  converge 
despite  the  singularity,  it  is  possible  to  improve  the  accuracy  of  some  of  the  "derived*  (to  be  defined 
later)  quantities  by  expanding  about  the  point  — /.  The  independent  variable  n  is  formed  by  the 
transformation: 

I 

fi  =  e^w  =  i  +  iy  (14)  ] 

where  /3,  a  positive  real  number,  is  defined  later.  Figure  2d  shows  the  mapping  from  2b. 


Fig  2  —  Domains  of  the  fluid  under  various  conformal  transformations:  a.  the  i  -  x  +  ty  plane;  b.  the  w  -  <t>  +  iii  plane;  c.  the 
mapping  of  the  fluid  into  the  interior  of  a  unit  circle;  d.  the  plane  appropriate  for  expansions  about  a  point  at  infinity.  The  point 
marked  $  marks  a  singularity  that  can  trouble  accurate  calculations.  It  is  infinitely  far  away  from  the  fluid  in  the  limit  of 
infinitesimal  waves;  it  approaches  the  free  surface  as  the  wave  amplitude  approaches  its  highest  value. 
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In  an  earlier  paper  making  use  of  this  transformation.  Wilting  (1975)  was  able  to  estimate  the 
location  and  some  of  the  properties  of  the  singularity  at  S.  For  the  wave  of  limiting  amplitude,  the 
highest  wave  in  water,  S  is  located  on  the  free  surface,  and  its  properties  are  known  (see  Lamb,  1932 
and  Grant,  1973). 


3.  TREATMENT  OF  THE  SINGULARITIES  AND  DEVELOPMENT 
OF  A  FOURIER  SERIES  SOLUTION 


A.  Basic  Relationships 


The  aim  here  is  to  find  a  solution  for  I)  (Eq.  (9))  that  satisfies  the  boundary  conditions  (Eqs. 
(10),  (11),  and  (13)).  We  split  the  solution  into  two  parts,  one  of  which  is  an  attempt  to  explicitly 
account  for  some  of  the  singular  behavior  at  S.  Thus, 


We  define  II 0  to  be 


ii  (?)  =  n0(O  +  «,<{) 

(15) 

ii0(0  =  00  +  ITq 

(16) 

!!,({>  -  »,  +  h,  . 

(17) 

,7  '  <  1  -  H 

3  TTa’ 

(18) 

where  X  is  a  constant  to  be  defined  presently.  Note  that  Eq.  (18)  satisfies  the  boundary  conditions 
(Eqs.  (10)  and  (ID),  but  not  necessarily  Eq.  (13).  On  the  unit  circle  Eq.  (18)  gives 


</0(<r) 


2X 


- r  ( 1  +  COS  u 

(1  +  X)2 


(19) 

(20) 


and  Eqs.  (19)  and  (20)  provide  the  relationship  between  X  and  <?o^) 

X  =  [l  -  </«3  (0)|  [l  +  Qo  (0)]. 


(21) 


In  the  numerical  work  X  is  to  be  specified,  sometimes  from  specifying  q0  and  using  Eq.  (21).  Hence, 
n0({)  is  determined  explicitly.  Note  that  for  the  highest  wave  the  singular  behavior  of  II 0  matches 
that  demanded  for  a  sharp  corner  along  a  free  surface. 

Specifically,  when  q0  vanishes,  X  is  unity,  and  the  singularity  in  Eq.  (18)  is  on  the  free  surface. 
Equation  (18)  goes  to: 


UQW  / 


ro(<J‘)  *  y  Injsi 


i 

sin  —  tr 


This  conforms  with  the  requirement  that  the  total  interior  angle  of  the  flow  is  120°,  for 

lim  0(<r)-  and  lim  # (<r )  —  -  . 

a~0*  0  <7 — tr  ® 


(22) 

(23) 
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As  g0  departs  from  zero  the  singularity  in  Eq.  (21)  moves  away  from  the  free  surface,  to  infinity 
when  <70  ”  1-  The  location  of  the  singularity  generally  conforms  to  the  estimated  locations  found  by 
Witting  (1975);  the  nature  of  the  singularity  conforms  only  for  the  highest  wave,  but  should  not  be  in 
too  much  error  for  almost  highest  waves.  In  any  case,  ft0  should  account  for  much  of  the  singular 
behavior  of  ft  at  S  for  high  waves. 

The  remaining  part  of  the  solution,  ft,,  is  represented  as  the  Fourier  series: 

ft,  s  «  L  £"•  (24) 

il  -  o 

The  boundary  condition  Eq.  (10)  demands  that  all  a„  be  real. 

The  components  of  ft,  on  the  free  surface  are: 


0,(<r)  = 

-  £  a„  sin  n<r 

(25) 

n  -  ! 

t,(<t)  = 

oo 

£  a„  cos  tut  . 

(26) 

n  —  0 

These  are  the  (real)  Fourier  expansions  of  9,  and  r,;  knowledge  of  either  as  a  function  of  o-  permits 
the  calculation  of  a„  through  various  techniques;  we  use  the  fast  Fourier  transform. 

Specifically,  Eq.  (13)  is  the  governing  equation  used  to  calculate  the  a„.  Only  the  interval  (0,jr) 
is  needed.  The  final  boundary  condition,  11,  is  that  q(-rr)  =  1.  Hence,  Eq.  (13)  has  the  extra  informa¬ 
tion  that 

1  =  ?3(0)  +  -  f  sin  0  sec  dtr.  (27) 

t tF1  2 

which  relates  F 2  to  q( 0)  and  fl(o-). 

Bernoulli’s  Law  gives  the  amplitude: 

a  -  y  FH\  -  q2m,  (28) 

and  by  quadrature  it  is  an  easy  matter  to  compute  a  profile  using  the  relationship 

I  pi»(<r)  _ 

dx  +  idy  = - r—~  sec  dv  (29) 

it  q(<r)  2 

that  arises  from  Eqs.  (8)  and  (9)  once  9  and  q  are  known. 

A  final  check  may  be  made  to  test  for  gross  programming  or  other  errors.  The  test  is  to  see 
whether  the  pressure,  initially  set  to  zero  by  Eq.  (7)  is  actually  zero  at  the  end.  This  (kinematic)  pres¬ 
sure  p  is: 

p-y(l-?2)-yF0'-l).  (30) 

Obviously,  that  p  defined  above  be  small  is  a  necessary  condition  for  an  accurate  solution;  however,  it  is 
not  a  sufficient  condition,  as  we  shall  see  later. 

B.  Derived  Relationships 

The  parameter  <?(0)  is  the  input  parameter  that  serves  to  distinguish  one  solitary  wave  from 
another.  We  call  those  other  parameters  basic  that  arise  most  directly  from  the  calculations  and  are 
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computed  to  the  highest  precision.  These  are  F 2,  a,  q,  and  9.  Derived  parameters  are  those  that  come 
less  directly  from  the  calculations,  or  cannot  be  computed  with  the  highest  precision.  The  first  of  these 
were  given  in  Eq.  (29),  namely  the  horizontal  and  vertical  coordinates,  from  which  a  profile  can  be  gen¬ 
erated. 

Other  derived  quantities  of  interest  are  defined  below.  The  solitary  wave  mass  M  is: 

M  =  J  -q  dx.  (31) 

where  tj  is  the  elevation  above  still  water  level.  A  parameter  closely  related  to  the  potential  energy  is 

f'v*.  (32) 

2  •s  oo 

Three  parameters  are  of  interest  at  the  flanks  of  the  solitary  wave.  An  expansion  in  powers  of  /j. 
(see  Eq.  (14))  leads  to  a  surface  profile  specified  by: 

x  =  <t>  -  x0  +  a^e^*  cos  p  +  0(e2fl*)  (33) 

t)  =  a sin  p  +  0(e2B*)  as  9  — ■  — oo.  (34) 

In  Eqs  (33)  and  (34)  p  satisfies 

F2  =  •  p  <  n/2.  (35) 

P 

The  three  parameters  of  interest  are  p,  x-0,  and  a,,  or  alternatively,  B,  defined  by. 

t)  —  Be  B'*'  as  l.xl  —  (36) 

The  horizontal  drift  2x0  is  identically  the  circulation  C  in  the  nondimensionalization  used  here. 
Knowledge  of  M  or  V,  and  C  is  sufficient  to  determine  other  integral  properties  of  interest,  using  the 
relations  given  by  Longuet-Higgins  (1974).  Once  profiles  have  been  computed,  it  is  easy  to  compute 
x0.  B,  and  either  Wor  V  or  both. 


4.  DEVELOPMENT  OF  NUMERICAL  SOLUTIONS 

A.  The  Solution  at  a  Fixed  Resolution 

Except  for  our  explicit  treatment  of  the  singularity  at  S  with  the  function  O0,  our  procedure  is 
similar  to  that  of  Yamada  (1957,  1958).  Table  1  gives  a  partial  description  of  the  steps  required  in 
forming  a  basic  iterative  solution  at  a  given  resolution,  i.e.,  value  of  N.  First,  a  value  for  </(0)  is 
specified.  This  identifies  a  unique  wave.  The  auxiliary  function  Jl0  is  chosen  so  that  q0( 0)  =  q( 0), 
which  determines  K.  A  total  of  N  +  1  Fourier  coefficients  a„  form  a  solution.  There  are  N  subdivi¬ 
sions  in  the  interval  0  <  rr  <  n.  Thus,  quantities  that  arc  functions  of  o-  are  determined  at  points 
(t„  =  nir/N  with  n  =  0,  1.  2 . N.  At  the  start  we  let  all  the  a„' s  be  zero. 

Steps  5  to  13  of  Table  I  outline  the  iterative  procedure  From  our  experience  the  most  delicate 
aspect  of  the  procedure  occurs  at  Step  7,  because  the  integral  appearing  in  Eq.  (27)  is  singular  at 
rr  =  n.  It  is  possible  to  show  that  ft  has  a  vertical  asymptote  at  <r  —  n  but  that  the  singularity  is  integr¬ 
ate  (we  find  no  reference  to  these  facts  in  earlier  work).  A  simple  approach  is  to  ignore  the  singularity 
and  evaluate  the  integrand  at  the  endpoint  by  using  the  limit 

lim  sin  9  sec  -^-<7  -  -  20’(ir)  (37) 

where  the  prime  denotes  the  first  derivative.  We  adopt  this  simple  approach.  Even  though  the  singu¬ 
larity  is  not  treated  elegantly,  the  error  should  be  small  for  N  sufficiently  large.  We  tried  to  develop 
methods  to  improve  our  treatment  of  this  singularity,  but  without  success. 
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Table  I  —  The  Numerical  Recipe 


Step 

Number 

— 

Step 

Equation(s)  and 
Methods  Used 

1 

Choose  i/(0),  N 

2 

Set  qn( 0)  =  q  (0) 

3 

Compute  (!„(£>  =  +  /ti,({) 

<! 

18,21 

4 

Set  a„  =  0  for  all  n  . 

This  completes  initialization 

5 

Compute  ) 

25 

6 

Compute  H  =  H „  +  w. 

7 

Obtain  F 2 

27,  Simpson's  Rule 

8 

Compute  a 

28 

9 

Test  to  see  whether  n  has  changed 
significantly  from  the  previous  iteration. 

If  not,  this  is  the  final  iteration. 

10 

Compute  q (it) 

13 

II 

Compute  7r(<r)  =*  In  q(ir)  -  rn((r) 

12 

Compute  a  new  set  of  a„'s 

26,  Fast  Fourier  Transform 

13 

Return  to  step  5  if  this  is  not  the  final 
iteration.  Continue  if  it  is  the  final  one. 

14 

Compute  ft,  q,  the  profile,  and 
other  derived  quantities 

Rather  than  specify  </(0),  one  can  just  as  well  specify  the  quantity  u>  defined  by 

w  =  I  -  fV(O).  (38) 

This  produces  some  slight  modifications  to  Table  1.  We  use  the  modified  procedure  to  produce  results 
that  can  be  compared  to  those  of  Longuet-Higgins  and  Fenton  0974)  and  Byatt-Smith  and  Longuet- 
Higgins  (1976) 

The  test  used  to  slop  the  iteration  process  is  that  successive  values  of  a  differ  by  less  than  10~6. 
Up  to  40  iterations  are  required.  An  example  of  the  convergence  process  is  shown  in  Table  2.  The 
results  display  features  found  in  every  case  examined:  (a)  The  higher  resolution  samples  (large  /V) 
start  much  closer  to  their  final  value  than  do  the  coarser  resolution  samples,  (b)  The  coarse  resolution 
samples  are  more  rapidly  convergent.  Fewer  iterations  are  usually  required  for  the  finer  resolution  cal¬ 
culations  The  slowness  of  convergence  for  /V  -  1440  in  Table  2  makes  it  risky  to  trust  the  last  digit. 
The  limiting  amplitudes  and  speeds  differ  significantly  between  one  resolution  and  the  other.  Why  this 
is  so  is  addressed  later. 
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Table  2  —  Convergence  of  Process 
to  Obtain  Solitary  Wave  Conditions  for  a 
Given  <u  The  value  of  <u  is  0.50. 


Iteration 

Number 

,V  = 

90 

.V  - 

1440 

a/h 

a / h 

F' 

1 

0  506102 

1  512203 

0480272 

1.460544 

2 

0.476858 

1  453715 

0  480275 

1  460549 

3 

0.465467 

1  430933 

0.480278 

1.460556 

4 

0.463457 

1.426915 

0480283 

1 .460566 

5 

0.465758 

1.43  151  7 

0.480290 

1.460579 

6 

0.469426 

1.438852 

0480297 

1.460595 

7 

0.472924 

1.445847 

0.480306 

1.460612 

8 

0.475618 

1.451236 

0.480314 

1 .460628 

9 

0.477404 

1.454809 

0.480321 

1.460642 

10 

0.478429 

1.456858 

0.480327 

1.460655 

II 

0.478917 

1.457833 

0.480332 

1.460664 

12 

0.479078 

1.458155 

0.480336 

1.460672 

13 

0.479070 

1.458139 

0.480339 

1.460677 

14 

0.478995 

1.457989 

0.480340 

1.460681 

15 

0.478909 

1.457818 

0.480342 

1 .460683 

16 

0.478839 

1.457677 

0.480342 

1 .460685 

17 

0.478790 

1.457580 

18 

0.478761 

1.457523 

19 

0.478748 

1.457495 

20 

0.478743 

1.457486 

21 

0.478743 

1.457486 

As  pointed  out  earlier,  a  necessary  condition  for  an  accurate  solution  is  that  the  pressure  found 
from  Eq.  (30)  nearly  vanishes.  Figure  3  displays  that  pressure,  also  for  tu  —  0.50.  The  computed  pres¬ 
sures  are  indeed  small,  never  exceeding  0.6  x  10-4,  even  at  the  coarse  resolution  N  *=  90.  They  are 
largest  near  u  —  180°,  the  flanks  of  the  solitary  wave.  A  striking  feature  in  Fig.  3  is  the  oscillation  of 
the  pressure,  with  an  alternating  sign  at  every  data  point  except  at  146°  to  148°.  This  oscillation  is 
characteristic  of  all  of  the  calculations,  except  for  isolated  points,  its  period  is  always  2  grid  points,  and 
so  is  related  to  the  calculations,  not  to  the  wave.  We  ascribe  the  oscillations  to  aliasing  in  the  fast 
Fourier  transform  algorithm,  though  it  is  possible  that  it  is  caused  by  some  other  numerical  artifact. 
Although  we  find  the  oscillations  unpleasant  to  look  at,  we  do  not  believe  that  they  indicate  a  serious 
problem. 

Figure  4  displays  the  behavior  of  the  computed  pressure.  On  a  log  scale  it  plots  the  envelope  of 
pressure  curves  like  that  of  Fig.  3,  for  N  =  180  and  N  =  1440.  The  curves  for  S  =  360°  and  for  N  = 
720  are  nested  within  the  pair  shown.  The  behavior  is  typical.  The  computed  pressure  increases  with 
<r,  is  always  small,  and  becomes  very  small  as  /V  becomes  large. 

B.  Extrapolation  to  Infinite  N 

Let  Q  be  the  exact  value  of  any  of  the  basic  or  derived  variables.  The  numerical  procedures 
described  in  Section  4A  estimate  Q  by  setting  some  finite  N  in  advance;  call  this  estimate 
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0\-  Presumably,  some  information  about  Q  should  be  present  in  a  sequence  of  Q ,v  that  is  absent  from 
any  individual  (>v.  We  guess  that  a  likely  behavior  is  that  of  a  power  law,  i.e., 

Qn  =  0  +  kNm.  <39) 

where  Q,  A,  and  m  are  unknown.  The  calculations  are  done  in  multiples  of  A ,  specifically,  A 
180,  360,  720,  and  1440.  If  Eq.  (39)  holds  exactly,  then  O  is  given  by: 

q  =  6/+1  Qj~\  ~  6/ 

Qj+i  +  Qj-i  -  2  Q} 

(see  Shanks,  1955).  It  turns  out  that  (Eq.  39)  appears  to  be  a  close  approximation,  in  that  the  esti¬ 
mates  of  Q  from  the  triads  360,  720,  1440  and  180,  360,  720  usually  are  much  closer  together  than 
C?i44o  is  to  either.  Moreover,  log-log  plots  of  {QN  —  Q )  versus  (TV),  which  must  be  straight  through 
the  points  A  =  1440,  720,  and  360,  also  pass  through  or  near  the  point  at  A  «  180  and  usually  near 
the  point  at  A  —  90.  We  believe  that  extrapolating  results  gives  a  closer  approximation  to  the  exact 
value  than  does  the  result  at  the  finest  resolution  (A  —  1440). 

C.  Loss  of  Accuracy  Near  the  Tails:  Connecting  the  Solution  to  an  Exponential  Decay 

The  use  of  the  fast-Fourier-transform  algorithm  demands  that  data  be  equally  spaced,  i.e.,  uni¬ 
form  intervals  in  <r.  As  we  have  seen,  this  does  not  produce  relatively  uniform  errors  (as  measured  by 
the  computed  pressures)  throughout  the  range  of  <r.  Neither  does  it  reach  very  far  out  on  the  flanks  of 
the  solitary  wave;  values  of  rjpjp  vary  between  0.044  (<o  =  0.1)  to  0.029  (to  =  1.0).  For  weak  waves, 
this  is  a  substantial  fraction  of  the  amplitude.  For  all  waves  we  match  the  Fourier  transform  solution  to 
an  exponential  falloff  at  the  tails  using  Eqs.  (33)  and  (34).  The  solutions  are  matched  at  <r  =  177°  for 
A  >  180. 

5.  SAMPLE  COMPARISONS  WITH  OTHER  WORK 

Longuel-Higgins  and  Fenton  (1974),  by  a  method  totally  different  from  ours  (an  expansion 
method),  use  &>  to  specify  wave  properties.  They  cite  5-place  accuracy  in  wave  speed  and  amplitude  up 
to  oi  =  0.75.  Through  this  region  (and  somewhat  beyond,  but  where  comparisons  can  be  made  only  to 
4-place  accuracy)  near-perfect  agreement  between  their  results  and  ours  exists.  Figure  1  shows  agree¬ 
ment,  but  does  not  indicate  how  remarkable  the  agreement  actually  is.  We  claim  about  5-place  accu¬ 
racy  from  <o  =  0.45  to  breaking.  Figure  5  shows  the  wave-speed  results  as  a  function  of  A  for  to  * 
0.45  and  u>  =  0.75  along  with  the  results  of  Longuet-Higgins  and  Fenton  (1974).  These  wave  strengths 
lie  at  the  ends  of  the  range  where  both  theories  claim  about  5-place  accuracy  (six  significant  figures  in 
F2). 

It  is  evident  from  the  figure  that  the  limiting  values  are  close— between  0  x  10-5  and  3  x  10~5.  It 
is  also  evident  that  the  extrapolation  procedure  helps  a  little  when  co  —  0.75  and  helps  greatly  when  co 
—  0.45.  This  is  characteristic.  The  extrapolated  value  of  all  parameters  lies  closer  to  the  A  —  1440 
value  the  closer  to  is  to  unity.  The  calculations  over  the  range  0<co^0.75  agree  with  those  of 
Longuet-Higgins  and  Fenton  (1974)  for  all  variables  given  by  them.  For  co  —  0.80  and  above  there  are 
small  (but  mathematically  significant)  differences  in  some  of  the  solitary  wave  parameters  (see  NRL 
Report  8505  for  details).  There  our  calculations  agree  more  closely  with  the  more  recent  work  of 
Byatt-Smith  and  Longuet-Higgins  (1976).  The  detailed  comparisons  with  high  solitary  waves  are 
deferred  to  NRL  Report  8505. 

Given  the  favorable  comparison  with  work  using  an  independent  method,  the  question  remains 
why  the  comparison  is  unfavorable  with  other  numerical  methods  (see  Fig.  I).  To  provide  some 
insight,  we  reduced  A  to  resemble  the  resolution  of  Yamada  (1957,  1958),  whose  numerical  method 
most  closely  resembles  our  own,  i.e.,  he  used  a  Fourier  Series  method  of  numerical  solution. 
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Fig.  5  —  Calculation  of  the  wave  speed  for  (a)  o>  -  0.75  and  for  (b)  w  -  0.45.  The 
points  show  the  speeds  resulting  from  calculations  at  the  various  resolutions 
specified  by  <V.  The  line  is  the  limiting  value  from  Eq.  (40)  from  estimates  at  JV  — 
360,  720  and  1440.  The  error  bars  show  the  speeds  found  by  Longuet- Higgins  and 
Fenton  (1974),  who  use  a  totally  independent  method.  They  give  F  to  five  decimal 
places;  the  error  bars  bracket  the  range  0.000005  from  their  stated  values. 


Specifically,  we  select  N  =  12,  24,  48,  90,  180,  360,  720,  and  1440  for  <?(0)  -  1/V3  (Yamada,  1958) 
and  for  <y(0)  =  0  (the  highest  wave;  Yamada,  1957).  Yamada  used  N  —  12,  probably  the  limit  for  the 
time.  Figure  6  draws  the  comparison.  Yamada’s  result  and  ours  for  /V  =  12  are  comparable  (his  grid 
system  and  numerical  procedure  differed  somewhat  from  ours,  however,  and  the  agreement  is  imper¬ 
fect).  We  conclude  from  the  figure  that  our  results  are  consistent  with  Yamada’s,  but  that  N  =  12  is 
insufficient  to  produce  even  2-place  accuracy  in  a  or  F2  —  I . 

Many  error  estimates  of  numerical  work  on  the  solitary  wave  problem  are  based  on  the  smallness 
of  the  pressure  as  an  indicator  of  errors  in,  say,  the  amplitude.  In  the  example  of  Table  2  and  Fig.  3,  w 
—  0.50,  the  pressure  residuals  are  small,  but  errors  in  the  speed  and  amplitude  are  not  so  small. 
Specifically,  for  N  *»  90  Table  2  indicates  an  error  of  0.0016  in  amplitude  and  0.0032  in  F2,  although 
the  residual  surface  pressures  never  exceed  0.00006  and  are  much  smaller  over  most  of  the  range  of  a 
(Fig.  3). 

Figure  7  shows  residual  pressures  and  estimated  errors  in  a  and  F1  for  u>  =  0.45.  We  see  that  the 
estimated  error  exceeds  the  surface  pressure  computed  at  all  parts  of  the  wave  by  two  orders  of  magni¬ 
tude  or  more;  the  excess  is  greater  than  three  orders  of  magnitude  for  a  central  point  (<r  =  90°).  In 
addition,  the  slope  of  the  estimated  error  differs  somewhat  from  that  of  the  computed  pressures  (which 
themselves  are  parallel).  This  means  that  a  simple  multiplicative  factor  is  insufficient  to  relate  com¬ 
puted  surface  pressure  to  errors  in  solitary  wave  properties.  We  find  that  residual  pressures  are  always 
much  less  than  estimated  errors  in  wave  speed  and  amplitude. 
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Fig  6  —  Calculation  of  the  wave  speed  and  amplitude  for  (a)  <?(0)  — 
1/V3  and  for  (b)  q( 0)  —  0,  the  highest  wave  The  open  points 
display  the  results  calculations  at  the  various  resolutions  specified  by 
N.  The  line  is  the  limiting  value  from  Eq.  (40)  from  estimates  at  /V 
*■  360,  720  and  1440  The  error  bars  show  the  data  from  (a) 
Yamada  (1958)  and  (b)  Yamada  (1957)  Calculations  with  resolu¬ 
tions  around  N  *  12  fail  to  give  solutions  accurate  to  two  decimal 
places  in  F2  or  a!  h 
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Fig.  7  —  Estimated  errors  for  a,  -  0.45.  The  uppermost  curve  estimates  the  error 
in  F1.  The  error  in  a/fi  is  one-half  of  this.  The  estimate  is  formed  by  taking  the 
data  at  N  -  350,  720,  and  1440  and  applying  E q.  (40).  This  ensures  that  the  last 
three  points  fall  on  a  straight  line;  that  the  first  two  also  fall  on  the  same  line  indi¬ 
cates  that  the  form  of  the  extrapolation  defined  by  Eq.  (40)  is  a  good  one.  The  data 
shown  with  open  symbols  give  the  residual  pressures  at  particular  points  on  the 
profile  of  the  solitary  wave.  It  is  evident  that  the  residual  pressures  are  poor  indica¬ 
tors  of  the  errors  in  speed  and  amplitude 


6.  SUMMARY  AND  CONCLUSIONS 

We  have  refined  the  numerical  method  developed  by  Yamada  (1957,  1958)  to  compute  the  pro¬ 
perties  of  solitary  waves.  The  major  refinements  are:  extending  the  resolution  of  the  calculations  to 
the  limit  of  our  computer,  probing  into  the  nature  of  the  behavior  of  some  of  the  resultant  parameters 
as  a  function  of  resolution,  and  extrapolating  to  the  limit  of  fine  resolution.  In  addition,  an  auxiliary 
function  is  included  to  mitigate  problems  caused  by  a  singularity  outside  the  fluid  near  the  crest. 
Unlike  some  more  recent  work,  we  went  to  the  Fourier  series  method  of  solution  for  two  reasons:  (a) 
the  fast  Fourier  transform  algorithm  permits  computation  involving  very  large  numbers  of  Fourier 
coefficients,  and  (b)  the  properties  of  Fourier  series  are  well  known  and  favorable.  Specifically,  the 
Fourier  series  is  convergent  even  for  functions  which  are  singular  at  a  point  on  the  unit  circle. 

The  agreement  with  calculations  using  nonnumerical  methods  is  impressive  Where  comparisons 
of  speed  or  amplitude  can  be  made  to  five  significant  figures,  there  is  disagreement  by  no  more  than  3 
x  10  \  except  for  very  high  waves.  There  is  no  reason  to  suspect  that  the  numerical  methods  used 
here  deteriorate  for  very  high  waves  (details  are  in  NRL  Report  8505).  This  is  the  first  numerical  treat¬ 
ment  of  solitary  waves  which  agrees  with  independent  methods  to  anywhere  near  the  four  decimal 
places  usually  stated. 
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What  went  wrong  with  the  other  numerical  methods?  In  the  case  of  Yamada  (1958)  we  have 
shown  that  the  number  of  Fourier  coefficients  used  is  insufficient  to  produce  highly  accurate  results. 
Moreover,  the  use  of  surface  pressure  as  an  indicator  of  errors  in  wave  amplitude  and  speed  can  be 
very  misleading— by  large  factors  (orders  of  magnitude).  We  suspect  that  the  same  problems  may 
occur  in  other  numerical  treatments,  i.e.,  insufficient  resolution  and  the  use  of  the  residual  surface 
pressure  to  provide  an  error  estimate.  In  addition,  the  singularity  at  o-  =  180°  may  cause  unrecog¬ 
nized,  perhaps-large  errors. 
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